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Abstract. Inspired by problems in biochemical kinetics, we study statistical properties of an 
over-damped Langevin process whose friction coefficient depends on the state of a similar, 
unobserved process. Integrating out the latter, we derive the long time behaviour of the 
mean square displacement. Anomalous diffusion is found. Since the diffusion exponent 
can not be predicted using a simple scaling argument, anomalous scaling appears as well. 
We also find that the coupling can lead to ergodic or non-ergodic behavior of the studied 
process. We compare our theoretical predictions with numerical simulations and find an 
excellent agreement. The findings caution against treating biochemical systems coupled 
with unobserved dynamical degrees of freedom by means of standard, diffusive Langevin 
descriptions. 
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1. Introduction 

Single molecule kinetics has come within reach of biophysical experiments 03 El El, and 
theoretical and computational tools for analysis of such processes have experienced a 
corresponding growth 0] [51 |6l Ul HO- h is clear that the combinatorially large number of 
microscopic steps involved in even the simplest of biochemical events makes their rigorous 
stochastic treatment difficult. For example, gene expression, often modeled as a single-step 
mRNA creation, in fact, includes transcription-factor-DNA binding, polymerase recruitment, 
transcriptional bubble formation, and multiple elongation steps, each of which is a complex 
process on its own. 

Therefore, any theoretical analysis of stochastic biochemical processes necessarily 
involves coarse-graining: identifying a small subset of dynamical variables that are modeled 
explicitly, while agglomerating the rest into an effective behaviour. Such coarse-grained 
dynamics is often modeled using the master equation or the Langevin approaches, which 
require Markoviness or white-noise random forcing. Both of these assumptions are, 
generally, flawed, and quantitative corrections have been worked out in certain cases [|9l[T0l. 
Less explored is the possibility that internal degrees of freedom can introduce qualitative 
differences, such as long-range temporal correlations among state transitions and non- 
diffusive behaviour of the observed quantities. 

A well-studied example shows that this is possible for a random walk on a discrete lattice. 
In Ref. [11], Weiss and Havlin analyzed a two-dimensional diffusion model, known as the 
comb model. There dynamics along the y coordinate is unlimited, while motion along x 
is allowed only when y = 0. This results in (x) = and (x 2 ) oc \/t, that is, in a sub- 
diffusive motion of x. This model is hardly realistic in a biochemical context due to the 
discontinuous dependence of the diffusion coefficient on y. However, it is plausible that 
diffusive dynamics of a real biological or chemical variable in the state space or in the 
physical space depends on unobserved, decimated variables in some other non-trivial way. 
For example, in a chemotaxing E. coli, the number of unobserved signaling proteins CheY- 
P is coupled to the distribution of times a bacterial motor rotates counterclockwise, and the 
bacterium swims straight. For a fixed concentration of CheY-P, obtained by modifying the 
chemotaxis network, the distribution is essentially exponential |fl"2l. resulting in a regular 
diffusive motion of the bacterium. But in a wild-type bacterium, as the number of CheY-P 
fluctuates (diffuses in the number space even for a constant external signal), the distribution 
becomes a power law, and bacteria exhibit super-diffusive real-space motion. While not true 
in this particular system, the distribution of clockwise rotation times could have been strongly 
coupled to CheY-P as well. This would have resulted in a power law distribution of times that 
the bacterium spends reorienting itself without moving forward, and hence in its sub-diffusive 
motion. In both cases, neglecting the CheY-P fluctuations and describing bacterial motion as 
a normal diffusion is qualitatively wrong. 

In this paper, we abstract out the detailed biology and explore these types of phenomena 
from the point of view of statistical physics. We derive the properties of a diffusion 
process, for which the diffusion coefficient depends on the state of another, unobserved, 
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variable. We show that, quite generally, such dependence leads to anomalous diffusion of the 
observed process, suggesting that traditional stochastic approaches may fail, and that more 
thought should be given to modeling stochastic phenomena in complex interacting systems, 
in particular in biophysics. 



2. The Model 



Our model is described by two variables x and y, which may represent, in particular, 
concentrations of two interacting chemical species. The x is considered as the observed 
quantity, while y is assumed hidden (i.e., unobserved). Particles of both species can be 
created and destroyed, which results in an over-damped diffusive motion of the system in the 
concentration space (we disregard the directional drift for simplicity, but it can be reintroduced 
easily). We assume that the diffusion of x is ^-dependent. That is, 

d y 1 (+\ m 
— = — 77(f), (1) 

at 7 y 

dx C (y) , . 

at j x 

Here 7 X , ^ y are the effective friction coefficients (assumed to be homogeneous) corresponding 
to x, y; i] (t) and £ (t) are independent, zero-mean white noise forces such that 

( v (t)r ] (t , )) = 2D yl 2 S(t-t / ), (3) 
(£(t)£(0> = 2D xl 2 J (t - 1') . ( 4 ) 



The idea of multiplicative noise was explored before. In [13] the authors considered the 
case of a Langevin process (not over-damped), in which a function of the velocity serves 
as a "filter" for the white noise. In |fT4j the multiplicative noise enters as a random friction 
coefficient; however, the distribution of the random friction is independent of time. Similarly 
in lfT5l the noise enters as a random mass with distribution which is again independent of time. 
Our model is different from the above mentioned models since the distribution of the random 
coupling parameter C (y) is time dependent; moreover it depends in an arbitrary way (through 
the function C(y)) on another variable, which diffuses and hence has long range temporal 
correlations. Further, the coupling considered in our model does not introduce a directional 
bias since it is multiplied by a white noise £, which takes both positive and negative values. 
The PDF of y is that of a normal diffusion 

1 (y-vo) 2 

p(y,t\y ,0)= == e «^ , (5) 

V^Dyt 

where the initial condition is y(t = 0) = yo. 

The dynamics of the mean square displacement (MSD) (x (t) 2 ), where (■ ■ ■) stands for 
the average over the white noises 77 and £, can be derived. Formally integrating Equation © 
and substituting it in the expression for the derivative of x 2 yields, 

dx(t) 2 2C{y) rt 



dt ^ x 



k{t) I C{y{t'))at')dt'. (6) 
^0 
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Figure 1. Typical trajectories for different forms of the coupling function C(y). We set 
D x = D y = j x = j y = A = 1. In the left panel we used the repressive coupling of Equation 
(0 with a = 0.625 to illustrate the suppression of diffusion in x. The central panel shows 
the case of a = in Equation ©, namely decoupled Langevin processes resulting in normal 
diffusion. The right panel shows the case of enhanced diffusion in x due to coupling of the 
form of Equation ( [ToT l with j3 = 0.75. In each of the panels the left inset shows the observable 
a; as a function of time and the right inset shows x and y scaled equally. All trajectories start 
at x = y = 0. The time step is set to 1 and the total duration of each trajectory is 10 6 . 



Averaging over the noise £(t) yields the dynamics of the MSD of x conditional on y(t), 

"<*y> = 2DMvm (7) 

To get the marginal expectation (x (t) 2 ), we now average the conditional expectation over y, 
which is distributed as in Equation ©: 



d(x(t) 2 ) 2D X f°° - C-»q>' , NN2 , 
at yjAxDyt J-oo 

The function C(y) may take different forms for different systems. The simplest case is 
when the dynamics of x is independent of y and C (y) = C = const. Substituting this in 
Equation © yields the expected trivial result (x 2 ) = 2D x C 2 t. 

Another scenario is that x can evolve in time only for a given range of y values 
(\y\ < yi), which resembles the discrete comb model of Weiss and Havlin IfTTTl . Notice, 
however, that the comb model is based on geometric constraints, i.e., the teeth, while in 
our case the coupling between the two process is not due to the topology, but rather to the 
physical nature of the processes. The similarity arises due the fact that in both models the 
first passage time distribution in the infinitely long y axis, dominate the dynamics. Indeed, 
substituting C(y(t)) = CO (y\ — y (£)) 6 (yi + y (t)), where C is a dimensionless constant, 
and Q(y) is the Heaviside Theta function, into Equation ©, we see that, for t y\f (12D y ), 
(x (t) 2 ) ~ \ft in agreement with IfTTTl . If C{y) falls off exponentially as y — > oo, the same 
sub-diffusion exponent is recovered. 

A more interesting case is when, at large y, C(y) falls off, but not too sharply. We 
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consider a power law form, namely 

= T7WP <9) 

where A is a constant with the units of inverse length, and a is a dimensionless parameter. 
This is of a form of a repressive, cooperative Hill kinetics, which describes many biochemical 
processes lfT6l . Typical diffusive trajectories with this C(y), A — 1, and a = 0,0.625 are 
shown in Figured] 

Assuming that the behaviour of C(y) for large y (i.e., C(y) ~ \y\~ a ) dominates the 
t — > oo dynamics of (x(t) 2 ), a simple scaling argument suggests that (x (t) 2 ) ~ t l ~ a . 
However, this is clearly wrong for large a, suggesting that (x (t) 2 ) must pick up an anomalous 
scaling due to the y — > properties of C(y). In what follows, we derive the long time 
behaviour of (x (t) 2 ) in a more rigorous way. 

Equation ([8]) with C(y) as in Equation ©, and considering the long time limit 
(VQ/y/ADyt < 1) gives 



dt 




dixit) 2 ) AD X f ajtd^i e" 

+ / 7 — , ; : scn2 d v > W 



[l+(V4^L4y) Q ] ' J-^ [I + {^DylAyY] 



In the first integral, we approximate the integrand as a constant for t — > oo, and, in the second 
integral, we neglect 1 compared to (^4D y t \Ay\) a , thus obtaining 



dt Ay/W~i (AA 2 D y y a ^ \2 AA 2 D y t 

where T (a, b) = J^°T a ~ 1 e~ T dr is the incomplete Gamma function. Integrating Equation (fTT)) 
over t results in the long time behaviour of (x(t) 2 ) 

Jxit) 2 )- J DiV / t + J D 2 t 1 "°, (12) 

where D lj2 are constants depending on the model parameters D x , D y , a and A. This implies 
that, for a < 1/2, the long time behaviour is dominated by (x(t) 2 ) ~ t 1-0 , as the scaling 
argument suggests. However, for a > 1/2, the scaling argument breaks and (x(t) 2 ) ~ \ft. 
Note that when the C (y) falls faster than 1 / y/y, the diffusion exponent is exactly the same as 
in the case in which the dynamics of x is limited to a finite range near y = 0. 

The case of a = \ is special and the integral of Equation ( fTOl can be calculated exactly, 
yielding 



^) 2 ) _ r54 f 1 
dt 45 4A 2 D„t 



o,o,i, |,f / Vt 



(13) 



where G denotes the Meijer G function IfTTTl . The leading order term of the MSD is then 



(x{t) 2 ) ~ \/tlnt. (14) 

The coupling function C(y) depends only on the variable y which undergoes normal 
diffusion. Thus we can derive the probability distribution of C(y). In order to simplify 
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Figure 2. Leading order diffusion exponent v defined by (x(t) 2 ) ~ t u for the coupled 
stochastic processes model with different couplings between the diffusion of x and y, measured 
by the exponent n. For the case of hindered diffusion of x, Equation (0, k = a, while for the 
enhanced diffusion, Equation ( TTol l. k = —(3. Numerical simulations (points) and theoretical 
predictions (line) agree for both scenarios. 



the notation, in what follows we use C to denote the coupling function without explicitly 
specifying its dependence on y. When C is given by Equation ©, its probability density is: 

_ (1-C) 2 / Q 

P (c) = — e iDytA * c i j (15) 

aC^" (1 - C)~ A^/^Ly 

where < C < 1. In the above derivation, we set = 0) = for simplicity. We see 
that the distribution of the noise introduced by C is time dependent. Note that the temporal 
correlations of C are also important for the dynamic properties of x. 

So far we have considered only situations in which the motion of x was slowed at large 
y, but we can also consider the opposite scenarios, when large y promotes diffusion in x, as 
in[[T2l: 

C (y(t)) = \Ayf, f3>0. (16) 

This now resembles the Hill activation kinetics for low concentration of the substrate 
molecules ifToll . Typical trajectories for this form of the coupling function with f3 = 0.75 
are shown in Figured! 
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For coupling of this form, the dynamics of the MSD(Equation ©) yields 

MA^fe-^)^^, (.7) 
at V 71 " J-oo 



which, in the long time limit, gives d(x(t) 2 )/dt ~ t 13 , and 

J^W)^t P+l - (18) 
The distribution of the coupling function C in this case is given by 

,3-1 C 2//3 

where C > 0. 

To confirm our analytical results, we performed numerical simulations for the different 
cases considered above. In Figure [2l we present a comparison of the diffusion exponent v 
(defined by (x(t) 2 ) ~ t u ) versus the coupling parameter k (for the sub-diffusion scenario, 
Equation ©, k = a, and for the super-diffusion scenario, Equation (fT6l) . k = —j3) between 
the simulations and the analytical results. The simulations where done according to Eqs. (Q][2]) 
with 7^ = 1, D XtV = 1 and dt = I. We averaged the results over 10 4 trajectories, each of 
duration 10 7 . . . W 8 dt. 

A simple linear regression to log (x(t) 2 ) = v\ogt + const was performed to estimate 
v. Since the standard parameter errors obtained for the regressions were negligible, the error 
bars of v were estimated from the variability of the fitted values as we changed the domain 
of t, for which the fits were performed. Figure [2] shows a clear agreement between our 
theoretical results and the simulations. Note that, in certain cases, the convergence to the 
leading behaviour of (x(t) 2 ) as t — > oo is slow since the difference between the exponents of 
the leading and the sub-leading terms is small. This slowness determined the lengths of the 
simulations. 



3. Time Averaged Mean Square Displacement 

There are many models of anomalous diffusion, including a time dependent friction coefficient 
in the Langevin equation [|T8|| , continuous time random walk (CTRW) |[T9ll , fractional 
Brownian dynamics 11201 . fractional Langevin dynamics ||2"T1 l22l l23l and Langevin dynamics 
with coloured noise |f24l|. to name a few. For a new model resulting in an anomalous diffusion, 
it is important to see if it can be reduced to one of these more familiar constructions. For 
example, the t — > oo behaviour of the original comb model is equivalent to CTRW ll25l with 
a power-law tail of the distribution of the times between successive jumps along x. However, 
in our model, the analogy is not as straightforward since the continuous dynamics of y induces 
temporal correlations among successive motions along x. 

To understand the relation of the coupled diffusion model to the others in the literature, 
we note that all of them yield the same behaviour for the ensemble averaged MSD in the 
long time limit. However, they still differ from each other in the short time behaviour, the 
shape of the distribution, and even in the long time behaviour of time averaged quantities 
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(for example, the CTRW exhibits ergodicity breaking 1126*1 ). In particular, the time averaged 
MSD (TAMSD) is an important property (it is the TAMSD that is observed in typical single 
molecule diffusion experiments in biological systems (31 |27l, and the number of recorded 
trajectories is often insufficient to estimate ensemble averages). 
The TAMSD is defined as 



which averages the squared displacement of a particle in time A (the time lag) over a single 
particle trajectory of duration t. For CTRW, the TAMSD is a random quantity and even 
its ensemble average still exhibits aging, that is, dependence on the measurement duration 
||28l |29l t in Equation (l"2"0l . On the contrary, for the fractional Brownian and Langevin 
dynamics, the TAMSD converges to the ensemble average MSD for long times [|30l . 

We investigated the behaviour of the TAMSD in our model with repressive coupling 
numerically. We find that, when the scaling argument holds, namely for a < 1/2 (see 
Equation ©), the TAMSD is not a random quantity, but it still shows aging, as we would 
expect for Langevin dynamics with a time dependent friction. On the other hand, when 
a > 1/2, and the diffusion exponent is v — 1/2, the TAMSD shows a similar behaviour 
to that of the CTRW (28). In Figure [31 we show the TAMSD versus the time lag A for 
a = 0.75 and a = 0.25. Each line shows the TAMSD for a single trajectory. The solid red 
lines represent trajectories of duration t = 10 5 , and the dashed blue lines are for t = 10 6 . 
All lines show a linear dependence on A, but with different coefficients. For a = 0.25, the 
TAMSD lines converge as the trajectory duration grows (the dashed blue lines essentially 
overlap), while for a = 0.75, the lines remain random. This is a clear indication of ergodicity 
breaking in our model for a > 0.5. Further, this analysis suggests that the coupled Langevin 
processes model stands as its own class among other anomalous diffusion models, exhibiting 
time-dependent diffusion coefficient Langevin dynamics for certain forms of coupling, and 
some aspects of ergodicity-breaking CTRW for others. 

4. Discussion 

In this article, we introduced a model of coupled diffusive processes, where the diffusion 
of the observed variable x is coupled to the value of a hidden variable y. We showed that 
the dynamics of x exhibits anomalous diffusion for every considered form of the coupling 
between the variables. Depending on the nature of the coupling, the motion of x can be sub- 
or super-diffusive (and even super-ballistic, as is the case of a frictionless particle subject 
to a white noise). Further, even for an arbitrary "strong" (such that the dynamics of x is 
limited to a small range of y values) repressive xy coupling, the x diffusion exponent v is 
limited from below by | (anomalous scaling), so that localization of x is impossible. Even 
though the long-time ensemble-averaged behaviour of our model is similar to that of many 
others describing anomalous diffusion, the model does not reduce to any of them, exhibiting 




(20) 
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Figure 3. The time averaged MSD S 2 (A, t) versus A for a = 0.25, 0.75 (left and right 
panels, respectively). All other parameters are the same as for the previous figures. We used 20 
trajectories of duration t = 10 5 (solid red lines) and 20 trajectories of duration t = 10 6 (dashed 
blue lines). The straight lines with slope 1 in the log-log scale reflect a linear dependence of 
the TAMSD on the time lag A. For a = 0.25, the TAMSD of different trajectories converges 
as the measurement time increases (the blue lines collapse, and the red do not), namely a time 
average of a single particle behaves like the ensemble average. However, for a = 0.75, there 
is no such convergence, just like for CTRW (the coefficient of proportionality between S 2 and 
A varies from trajectory to trajectory). In both cases, the ensemble average of the TAMSD 
decreases as the measurement time increases, indicating aging. 

an effective time-dependent diffusion coefficient, aging, and ergodicity breaking for different 
values of its parameters. 

The anomalous scaling and the ergodicity breaking appear for "strong" xy coupling (i.e., 
large a). This is because, for a < 1/2, motion of particles away from y = contributes 
substantially to the ensemble averaged MSD of x. On the contrary, for a > 1/2, only motion 
near y — is important, namely the first passage time in an infinite one dimensional system 
(the diffusion process in y) plays a major role in the dynamics of the observed quantity, x. 
It was also verified numerically that, for a > 1/2, where the dynamics of x is dominated 
by motion in a narrow range near y = 0, the propagator is non-Gaussian as expected from a 
CTRW-like renewal process. A similar result holds for the super-diffusive regime, j3 > 0. On 
the other hand when the motion at any y is important (0 < a < 1/2), the process in x has 
long-time correlations and the propagator takes a Gaussian form (similar to that of a diffusion 
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process with a time dependent diffusion coefficient). 

While important in its own right, the coupled diffusion model raises its most important 
questions in the biological domain. Unobserved dynamical quantities lead to anomalous 
diffusion in E. coli chemotaxis lfT2l . or in mRNA diffusion in cells 0. Further, some of 
the best established models of cellular regulation involve coarse-graining of dynamics. For 
example, in some models of the E. coli lac operon, the lac repressor itself is an unobserved 
variable OTI . which is coupled to the speed of production of the lactose permease and the 
lactose-utilizing enzyme and, through them, to the import and the degradation of lactose in 
the cell. Since any coupling may lead to anomalous diffusion and in some cases even to 
ergodicity breaking, it begs the question whether relying on the common Langevin or master 
equation analysis of stochasticity of the lac repressor or other regulatory circuits, such as the 
A-phage [@J[8]|, 01, and others, is rigorously justifiable. 
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